/***************************************************************************
 *   Copyright (C) 2008 by Ivan Blagoev Topolsky   *
 *   ivan.topolsky@genebio.com   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU Library General Public License as       *
 *   published by the Free Software Foundation; either version 3 of the    *
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this program; if not, write to the                 *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream>
#include <iomanip>

// BOOST
#include <boost/lexical_cast.hpp>

// Captain-bol
#include <spectrumpeakrefdualptr.h>
#include <chargemapper.h>

// OLE/COM support
#include "comutil.h"
#include "disphelper.h"

#include "analyst.h"


#ifndef DEBUG_DUMP_CSV
#define DEBUG_DUMP_CSV 0
#endif


namespace xeno {
namespace freud { // JOKE freud's an analyst.

using namespace std;	// STL
using namespace boost;
using namespace xeno;
using namespace bol;	// Captain bol

// FIXME the following isn't portable between bytesex. But as currently OLE/COM only runs on bigendian x86 architecture it doesn't matter much (minus DarWine on PPC).
typedef union {
	long long refid;
	struct {
	// BUG as we have to cram everything inside a 64bit REFID I currently use 16bit shorts - BUT the actual numbers returned by analyst are 32bit longs
	// overflow could happen, BUT I have yet to find a WIFF file with more than 64k sample or 64kperiods
	// Given that Analyst loads whole periods into RAM, I don't even think that 64k experiments or 64k cycles are imaginable
	// FIXME find some idea to fix this. Maybe by multiplying everything by "num..." and thus having a continuous numbering scheme instead of a bit field
	// FIXME refid = sample + numsam*(period + numpe*(experiment + numxp*(cycle)))
		unsigned long cycle;
//		unsigned short cycle;
		unsigned short experiment;
		unsigned short period;
//		unsigned short sample;
	} th;
} REFID;

long long analyst::mkrefid() const {
	REFID r;
	r.th.cycle = curcy;
	r.th.experiment = curxp;
	r.th.period = curpe;
//	r.th.sample = cursam;

	return r.refid;
}
	
/** 
Hierachical structure of Wiff file :
- A Wiff file is a collection of data about a series of samples.
WARNING: The samples are counted from 1 to N (*not* 0 to N-1)

- A sample can contain several periods - different groups of experiment that were done at different times.
basically, that means it's possible to do completely different runs on the same sample
but normally only 1 period is used in most lab workflows.
WARNING: Whenever looking inside an experiment, Analyst freeze for a certain amount of time processing data
TODO ask PAB and see with the people of chemical analysis how the elution-time specific detection will be implemented

- for each period we have several "experiments" : specific settings at which data was collected
usually most lab workflow use "IDA" :
- during the whole elution of the sample through the LC column, the spectrometers constantly does a simple analysis of the elute
(using plain MS or MRM) and continuously store the data.
this data and the corresponding settings are grouped as "experiment #0"
- whenever the "experiment 1" picks up something interesting, more detailled analysis will be performed on the elute (some king of MS/MS). 
Usually these analisys are performed using several different conditions (varying the collision energy).
each set of data for a given condition is stored in 1 experiment. 
(most lab workflow will have 3 different collision energies, stored as "experiment #1", "experiment #2" and "experiment #3")

- the repeated mesures of an experiment are called "cycles". There are as many cycles as there are spectra per experiment or points per chromatogram
*/


/*
 * Data types used by Analyst
 */

// scantype
typedef enum scan_enum {
	// from Analyst CookBook
	SCAN_Q1 = 0, QCAN1_MI, SCAN_Q3, SCAN_Q3_MI, SCAN_MRM, SCAN_PREC, SCAN_PROD_ION, SCAN_NEUT_LOSS, 
	// hapilly stolen^H^H^H^H^H^H borrowed from mzWiff
	SCAN_TOF_MS, SCAN_TOF_MS2, SCAN_TOF_PREC, SCAN_EPI, SCAN_ER, SCAN_MS3, SCAN_TDF, SCAN_EMC, SCAN_EMS, SCAN_ELSE, SCAN_UNKNOWN
} SCANENUM;

static const struct {
	const char *name;
	const int level;
} ScanType[] = { 
	// from Analyst CookBook
	{"Q1 Scan",1},
	{"Q1 MI",1},
	{"Q3 Scan",2},
	{"Q3 MI",2},
	{"MRM",2},
	{"Precursor Scan",1},
	{"Product Ion Scan",2},
	{"Neutral Loss Scan",1},
	// shamelessly copied from mzWiff
	{"TOF MS",1},
	{"TOF MS2",2},
	{"TOF Precursor Ion Scan",1},
	{"EPI",2},
	{"ER",1},
	{"MS3",3},
	{"TDF",1},
	{"EMC",1},
	{"EMS",1},
	{"Else",1},
	{"Unknown",1},
	{NULL, 0}};

/**
 * this function takes care of the first steps of handling raw files :
 * - getting the necessary OLE/COM object
 * - setting up a RunLcmsms object to hold the results
 * 
 * @param ifname Filename to use
 */

void analyst::initfile(const char *ifname) {
	// initialise objects
	try {
		// Analyst.FMANWiffFileControl.1
		dhCheck( dhCreateObject(L"Analyst.FMANWiffFile.1", NULL, &wiff) );
		dhCheck( dhCreateObject(L"Analyst.FMANChromData.1", NULL, &chrom) );
		dhCheck( dhCreateObject(L"Analyst.FMANSpecData.1", NULL, &spec) );
	} catch (string errstr) {
		cerr << "Fatal error - can't create the main Analyst objects" << endl
			/// This software needs a couple of COM/OLE classes registered by a few DLLs
			/// As said in the help, you can either install them 
			/// Or manually register them
			<< endl << "Is the library correctly registered ?" << endl
			<< "- either run Analyst's installer" << endl
			<< "- or manually run :" << endl
			<< "\t> regsvr32 {drive}:\\{path}\\ExploreDataObjects.dll" << endl
			/**
			module "ExploreDataObjects" - objects for manipulation of file and data within
			contains classes :
			- FMANWiffFileControl	(open wiff files)
			- FMANWiffFile	(the wiff files themselves)
			- FMANChromData
			- FMANSpecData
			*/
			<< "\t> regsvr32 {drive}:\\{path}\\MSMethodSvr.dll" << endl
			/**
			module "MSMethodSvrLib" - parameters of experiment
			contains class :
			- Experiment
			 */
			<< "\t> regsvr32 {drive}:\\{path}\\ExploreDir.dll" << endl
			/**
			module "ExploreDirLib" - used to process data
			contains class :
			- CentroidOptions
			- ExploreProcessing
			 */
			<< endl << "details:" << endl << errstr << endl;
		exit(2);
	}

	/**
	NOTE Normally, according to the Analyst Cookbook, we open a file using the FileControl with something like this (Visual Basic :)
	
	DIM filecontrol AS FMANWiffFileControl = new FMANWiffFileControl
	DIM wiffile AS IFMANWiffFile = filecontrol.GetWiffFileObject(fname, 1)
	
	*BUT* surprisingly, lots examples in theCookBook use directly FMAN*Data objects,
	using  "->WiffFileName" property to get/set the filename
	and "->GetWiffFileObject" to get the wiff file object.

	mzWiff uses the same approach
	
	NOTE This shows that Analyst has probably an ugly internal design where property getter/setter aren't pure functions, but instead hook calls to do other behind the scene work. The ugly design could explain a lot of the weird behaviour (like Analyst not working over network).
	*/
	
	// opening 
	cout << "reading from : " << ifname << endl;
	try {
		dhCheck( dhPutValue(spec, L"WiffFileName %s", ifname) );
		dhCheck( dhPutValue(chrom, L"WiffFileName %s", ifname) );
		dhCheck( dhGetValue(L"%o", &wiff, spec, L"GetWiffFileObject") );
	}catch (string errstr) {
		cerr << "Can't openfile " << ifname << endl

			<< endl << "Are you trying to open a file from a network location ?" << endl
			<< "the Analyst API is known (according to the Cookbook) to have difficulties when openning Wiff file from the network." << endl
			<< "- either try to copy the file to a local drive and run again" << endl
			<< "- or first load the file into Analyst 1.4.x it self and then try again" << endl
			// Yes, the Analyst API is that bad

			<< endl << "Also, if running from a POSIX shell (like Cygwin's bash) :" << endl
			<< "- please use MS-DOS like paths with backslashes ( \\{path}\\{file}.wiff )" << endl
			<< "- instead of POSIX paths with slashes ( /{path}/{file}.wiff )" << endl
			// Analyst wasn't probably written entirely in C++. That would explain the problem with paths
			// It does partially work, though. But on a random fashion.
			<< "- and do NOT use extended POSIX-like features provided by Cygwin ( no symlinks, no mount points, etc. )" << endl
				
			<< endl << "details:" << endl << errstr << endl;
		exit(2);
	}


	cout << endl;
	// Normally, run flow continues with the next function
	// they are currently split to allow separate calling 
}


/**
 * normally this step comes right after the first
 * - it will display a list of all samples contained in a wiff file
 * - it will store the names in a string vector (samnames) for further usage
 */
void analyst::samplelist() {
	/// In a WIFF file, there can be more than one sample

	//
	// Show all samples
	//
	numsam = 0;	
	cout << "samples : ";
	try {
		dhCheck( dhGetValue(L"%d", &numsam, wiff, L"GetActualNumberOfSamples") );
		
		samnames.resize(numsam+1, "");
		cout << numsam << endl;
		
		for (int i = 1; i <= numsam; i++) {
			char *name = NULL;
			dhCheck( dhGetValue(L"%s", &name, wiff, L"GetSampleName %d", i) );
			
			samnames[i] = name;
			cout << " " << i << ":\t" << name << endl;
			
			dhFreeString(name);
		}
	} catch (string errstr) {
		cerr << "Error getting samples" << endl
			<< endl << "details:" << endl << errstr << endl;
		exit(2);	// TODO maybe find a more elegant solution when failing
	}
	cout << endl;
}


/**
 * dosample processes 1 single sample
 * @param sample	sample number to process (warning: samples are 1 to N, not 0 to N-1)
 * @param ptr	Write to a specific RunLcmsms instance. Other wise, if NULL, we allocate one.
 *
 * Also uses the following members:
 * wiff 	pointer to an instance of coclass FMANWiffFile
 * chrom	pointer to an instance of coclass FMANChromData
 * spec 	pointer to an instance of coclass FMANSpecData
 */

void analyst::dosample(const int sample, RunLcmsms *ptr, const bool skipms1, const bool skipms2, const bool dropms1, const bool dropms2, ExtractList* reflist) {
	// allocate an instance if needed...
		// *** EXPORT CAPTAINBOL
	run = ptr ? ptr : new RunLcmsms();

	/**
	For a given sample there is 1 or more "periods" :
	different time period when measure were done on the elute comming out of the LC column

	this make it possible to only record data during a couple of key moments 
	(for example at the precise time when some substances are known to elute)
	by using a couple of periods, rather having the whole MS reading running during the whole elution time.
	
	but currently only 1 period is used in most lab workflows.
	
	TODO ask PAB and see with the people of chemical analysis how the elution-time specific detection will be implemented
	
	There are function called :
	- 'GetActualPeriodAndCycleFromTime'
	- 'GetActualTimeFromPeriodAndCycle'
	which helps converting between absolute time and periods/cycles
	*/
	
	assert(sample <= numsam);
	assert(sample < samnames.size());
	
	cursam = sample;
	
	samnames[0] = samnames[sample];
	string& samname = samnames[0];
	
	initextract_map(reflist);
	
	try {
		dhCheck( dhGetValue(L"%d", &numpe, wiff, L"GetActualNumberOfPeriods %d", sample) );
	} catch (string errstr) {
		cerr << "Error getting samples" << endl
				<< endl << "details:" << endl << errstr << endl;
		return;	// TODO maybe find a more elegant solution when failing
	}

	// sample name
	cout << "sample " << sample << " - " << samname << " - Per:" << numpe;
	// *** EXPORT CAPTAINBOL
	run->setDescription(samname);
	for (int pe = 0; pe < numpe; ++pe) {
		curpe = pe; // for status, etc.
		numxp = numcy = 0;

		// Get number of experiments
		try {
			/// NOTE Analyst usually freezes when accessing the content of a period for the first time
			dhCheck( dhGetValue(L"%d", &numxp, wiff, L"GetNumberOfExperiments %d %d", sample, pe) ); // <- exactly here
			// TODO provide a nice status window during this time. Or at least a hour glass or something
		} catch (string errstr) {
			cerr << "Error getting periods" << endl
				<< endl << "details:" << endl << errstr << endl
				<< endl << "trying to continue with next period" << endl;
			continue;
		}

		if (numpe > 1) {
			if (pe == 0) cout << endl;
			cout << "Period:" << pe;
		}
		cout	<< " - Exp:" << numxp;


		// Get number of cycles
		try {
			/// NOTE Analyst usually freezes when accessing the content of a period for the first time
			dhCheck( dhGetValue(L"%d", &numcy, wiff, L"GetActualNumberOfCycles %d %d", sample, pe) ); // <- exactly here
			// TODO provide a nice status window during this time. Or at least a hour glass or something
		} catch (string errstr) {
			cerr << "Error getting periods" << endl
				<< endl << "details:" << endl << errstr << endl
				<< endl << "trying to continue with next period" << endl;
			continue;
		}
		
		cout	<< " - Cyc:" << numcy << endl;


		for (int xp = 0; xp < numxp; ++xp) {
			curxp = xp; // for status, etc.
		/** as it is an application for controlling an industrial spectrometer doing experiment, interesting data is split across two different place
		- settings, parameters, etc. In short : everything the scientist entered in the software before the machine started to work
		is stored into various structure of the MSMethodSrvLib DLL. 
		for example the collision energy for a given experiment (= series of data, see above) will be stored here.
		- all the actual readings that came from the detection unit of the mass spectrometer, 
		are acessible through the various structure provided by ExploreDataObjects
		*/
			SCANENUM scantype = SCAN_UNKNOWN;
			bool negpol = false;
			string desc ("");
			
			try {
				char *ptr = NULL;
#if 0
				/*
				For some reason "GetExperimentObject" fails to get the *actual* experiment
				And instead creates a new empty one (called "New experiment" and by default of + Q1Scan type)
				*/
				CDispPtr exp;
				
				dhCheck( dhGetValue(L"%o", &exp, wiff, L"GetExperimentObject %d,%d,%d", sample,pe,xp) );

				dhCheck( dhGetValue(L"%d", &scantype,	exp, L"ScanType") ); 
				dhCheck( dhGetValue(L"%b", &negpol,	exp, L"Polarity") );
				dhCheck( dhGetValue(L"%s", &ptr,	exp, L"Description") );
				
				// other interesting properties are : ScanMode 
#else
				/*
				Contrived way to get the same information using FMANSpecData's GetScanInfo method
				*/
				double _sec = 0.;
				dhCheck( dhCallMethod(chrom, L"SetToTIC %d,%d,%d", sample,pe,xp) );// look inside the TIC
				dhCheck( dhGetValue(L"%s", &ptr, chrom, L"DataTitle") ); // HACK this is the closest to a description I can get
				
				// dhCheck( dhGetValue(L"%e", &_sec, chrom, L"GetXValueInSec %d", 1) ); // get time of fist spectrum
				dhCheck( dhGetValue(L"%e", &_sec, wiff, L"GetActualTimeFromPeriodAndCycle %d,%d,%d,%d",sample,pe,xp,0) ); // get time of fist spectrum

				short _gettype = SCAN_UNKNOWN; // as we use call by ref (pointers) we need to have the exact type width
				short _getmode = 0;
				long _negative = 0;
				short _dumpq1 = 0;
				short _dumpq3 = 0;
				dhCheck( dhCallMethod(spec, L"SetSpectrum %d,%d,%d,%e,%e", sample,pe,xp,_sec,_sec) ); // get the first spectrum
				dhCheck( dhCallMethod(spec, L"GetScanInfo %&hd,%&hd,%&ld,%&hd,%&hd", &_gettype, &_getmode, &_negative, &_dumpq1, &_dumpq3) );
				scantype = (SCANENUM)(int)_gettype;
				// scanmode = _getmode;
				negpol = _negative;
#endif
				desc = ptr;
				dhFreeString(ptr);
			} catch (string errstr) {
				cerr << "Error getting informations about experiment" << endl
					<< endl << "details:" << endl << errstr << endl
					<< endl << "trying to continue with next experiment" << endl;
				continue;
			}

			cout << "\ttype:" << scantype 
				<< "\t" << (negpol ? "-" : "+")
				<< " " << ((scantype < SCAN_UNKNOWN) ? ScanType[scantype].name : ScanType[SCAN_UNKNOWN].name) 
				<< "\t\"" << desc
				<< "\"\tlevel:" << ((scantype < SCAN_UNKNOWN) ? ScanType[scantype].level : ScanType[SCAN_UNKNOWN].level);
			cout << endl;

			/** IDA
			as explained above LC MS/MS runs are usually setup with a mode that Analyst calls IDA :
			several "experiments" i.e.: several different type of scan are done during the elution.
		
			"experiment #0" - the first type of scan, is systematically done on all cycles
			it is used to quickly identify potential candidate for analysis in the elute.
			often it a simple MS^1 scan - which is simply used to detect precusors and their masses
			it can also be y MRM (see bellow) - with is used to recognise potentially known substance
		
			"experiment #1" and the following - are done only on interesting candidates determined by the previous one.
			they are MS^2 and are used to identify ions found on an MS^1
			or to analyse and confirm identity of potential candidates found with MS^2
			*/
			if (scantype == SCAN_MRM) {
				cout << "\tprocessing MRM..." << endl;
				if (xp > 0) 	cout << "\tshould be the first ?!?";
				
				if (! dropms1) {
					dospectra(&analyst::do_single_spectrum_mrm, sample,pe,xp, skipms1);
					cout << endl;
				} else
					cout << "\tuser skipping all" << endl;
				continue;
			}
			
			switch(ScanType[scantype].level) {
				case 1:
					// MS level 1 - usually the experiment #0 in an IDA
					// used to detect precursors
					if (xp > 0) 	cout << "\tthis scan should be the first ?!?";
					cout << "\tgetting precursor..." << endl;
					
					if (! dropms1) {
						dospectra(&analyst::do_single_spectrum_ms1, sample,pe,xp, skipms1);
						cout << endl;
					} else
						cout << "\tuser skipping all" << endl;
					continue;
				
				case 2:
					// experiment #1 and following : detailed analysis
					if (xp == 0) cout << "\tthere should be a level 1 scan here ?!?";
					cout << "\tgetting spectra..." << endl;
				
					if (! dropms2) {
						dospectra(&analyst::do_single_spectrum_ms2, sample,pe,xp, skipms2);
						cout << endl;
					} else
						cout << "\tuser skipping all" << endl;
					continue;

				case 3:
					// experiment #1 and following : detailed analysis
					if (xp == 0) cout << "\tthere should be a level 1 scan here ?!?";
					cout << "\tgetting MS3..." << endl;	
					// TODO implement MS^3 in RunLcmsms
				default:
					cout << "\t\tNot implemented yet !!!" << endl;	
					continue;
			}
		} // experiment
	} // period
}



/**
 * This step is called to process all spectra in all cycles for a given experiment
 *
 * Parameter defining Action :
 * @param effector decide wether we'll be processing parent MS or fragment MS^2
 * its a pointer to the corresponding function whitch will take care of the spectrum itself
 *
 * Parameters how to process the data
 * @param docentroid	FIXME not implemented yet
 * @param skip	do we process the data ? or do we skip it
 *
 * Global data used :
 * cur... : which experiment we are scaning
 */
void analyst::dospectra(const EFFECTOR_FUNC_PTR effector, const int sam, const int pe, const int xp, const bool skip, const bool docentroid) {
	/**
	- ChromData object gives a bird's eye view of everything that was hapenning during the whole elution time (chromatogram)
	- SpecData gives detailed data about what was detected at a specific point in time (spectra)
	- At each point in time in a chromatogram, there are spectra aquired
	- These spectra are addressed using *single precision* floats. When using doubles the "time" may end up pointing between spectra
	- we can pick the sum of spectra for a whole time period, but usually we will only want the data at the precise moments of the chromatogram

	WARNING data points are numbered from 1 upward (like sample)
	
	NOTE it's also possible to find spectra using GetActualTimeFromPeriodAndCycle and not use the chromatogram. But in that case we lose the check for empty spectra (shows as 0 on a TIC)
	*/
	int chpts = 0;
	try {
		dhCheck( dhCallMethod(chrom, L"SetToTIC %d,%d,%d", cursam,curpe,curxp) ); // get the data : overall chromatographe
#if DEBUG_DUMP_CSV
		string debugout ("out_" + boost::lexical_cast<string> (curxp) + "_chrom.csv");
		dhCallMethod(chrom, L"SaveToTextFile %s", debugout.c_str());
#endif
		
		dhCheck( dhGetValue(L"%d", &chpts, chrom, L"GetNumberOfDataPoints") );
		cout << "\tNumber of data points: " << chpts;

		double min = 0., max = 0.;
		
		dhCheck( dhCallMethod(chrom, L"GetXValueRange %&le,%&le", &min, &max) );
		cout << "\tThere's data between " << min << " and " << max << endl;

		if (numcy != chpts)
			cerr << "weird mismatch of cycles and chromatogram !!!" << endl;
		
		// scan the whole chromatogram -> the whole duration of the elution
		for (int i = 1; i <= chpts; ++i) {	// Warning: DataPoints are 1 to N (not 0 to N-1 like cycles)
			curcy = i-1; // Warning: cycles are 0 to N-1 (not 1 to N like DataPoints)
		// this loop parse the mass readings
		/** semantic consideration about MS measurment
		- inside a chromatogram, it's simply called a "data point" - Beware: 1 to N (*not* 0 to N-1) this contradicts "cycles"
		- inside the "msmethod" classes, it's called a cycle - Beware : 0 to N-1 (*not* 1 to N)
		- spectra aren't addressed directly but by time-frame - Beware : this is *always* specified in seconds (whatever the default time scale of the chromatogram is)
		*/
			double tic;
			dhCheck( dhGetValue(L"%e", &tic, chrom, L"GetDataPointYValue %d", i) );  // total intensity at this point in time / this cycle
	
			// are there at least some detected ions ?
			if (! (tic > 0.)) continue;	// <- that's why the chromatography codepath is faster !
			
			// (So we'll use the time reported by the chromatogram to get corresponding spectra)
			double t;
			dhCheck( dhGetValue(L"%e", &t, chrom, L"GetXValueInSec %d", i) );  // total intensity at this point in time / this cycle

			(this->*effector) (sam,pe,xp, t,tic, skip,docentroid);
		}
	} catch (string errstr) {
		cerr << "Error getting chromatogram" << endl
				<< endl << "details:" << endl << errstr << endl
				<< endl << "skipping experiment" << endl;
		return;
	}
}



/**
 * This step decodes 1 single parent spectrum
 *
 * Spectrum coordinates (FIXME this duplicates  the content of "cur???" globals) :
 * @param sam	current sample (also in "cursam" global)
 * @param pe	current period (also in "curpe" global)
 * @param xp	current experience (also in "curxp" global)
 *
 * Parameters about spectrum :
 * @param tic	tic of the spectrum (if known)
 * @param t	retention time of the  spectrum
 * 
 * Parameters how to process the data
 * @param skip	do we process the data ? or do we skip it
 * @param docentroid	ignored : MRM aren't continuous, they are just a collection of separate reactions monitoring
 */
void analyst::do_single_spectrum_mrm(const int sam, const int pe, const int xp, const float t, const double tic, const bool skip, bool docentroid) {
	/** MRM - Multiple reaction monitoring

	mostly like Reaction Monitoring (a single peak MS^2) which is made to quickly identify candidates that are potentially known substances :
	- based on their mass (q1 with fixed mass)
	- and on a specific fragment (q3 with fixed mass)
	
	an MRM is a collection of such Reaction monitoring, designed to look for multiple such candidates
	- it could be considered as a collection of multiple single-peak MS^2
		(that would be the most logical for bioinformatics processing)
	- or as a weird spectrum whose X axis aren't simply M/Z, but M/Z pairs (for q1 and q3)
		(that's what Analyst does. It has the advantage of being more compact to display on screen)
	
	This type of scan is useful when screening a sample for known substances.
	
	HACK currently we consider that an MRM shows precursors
	so the MRM is processed to create a "virtual" MS^1 spectrum of precursors
	it's not fundamentally wrong, as traditionnally those are represented as chromatograms
	(per each precursors, a collection of intensity over time)
	wheras we simply have the data packaged the other around
	(per each time slot, a collection of intensity over precursors' m/z)

	TODO at some time later once the mecanism is correctly added in CaptainBOL, we should use a special class
	PeakIntensityRM which can also hold additionnal informations about the MRM/SRM
	*/

	cout << "\t\tProcessing MRM " << curcy << "\tat " << t << "\ttic:" << tic;
	int sppts = 0;
	int scans = 0;
	string title("");
	
	try {
		char *ptr = NULL;
		dhCheck( dhCallMethod(spec, L"SetSpectrum %d,%d,%d,%e,%e", sam,pe,xp,t,t) ); // at the same exact time as the point on the chromatogram
		// note : whatever the default time scale in the wiff file, the "SetSpectrum" method *needs* seconds.
		
		dhCheck( dhGetValue(L"%d", &sppts, spec, L"GetNumberOfDataPoints") );
		dhCheck( dhGetValue(L"%d", &scans, spec, L"GetNumberOfScans") );
		dhCheck( dhGetValue(L"%s", &ptr, spec, L"DataTitle") );
		title = ptr;
		dhFreeString(ptr);
	} catch (string errstr) {
		cerr << "Error getting spectrum" << endl
			<< endl << "details:" << endl << errstr << endl
			<< endl << "skipping spectrum" << endl;
		return;
	}
	
	cout	<< "\tnum of points: " << sppts
//		<< "\tspectra:" << scans	
			// Well we just have to pay attention that no more than 1 scan are used in the current spectra
		 	// (unless combining several scan in one on purpose)
			// TODO optionally we could add a few asserts but it's not required
		<< endl;
	
	if (scans > 1) {
		cerr << "we should use data obtain from only 1 spectrum !" << endl;
	}
	
	// NOTE no centroiding as MRM aren't continuous

	/*
	* HACK sum together intensities for precursors with similar mass
	* TODO once PeakIntensityRM are implemented this should be deprecated and we simply build an instance of PeakIntensityRM using Q1 and Q3
	*/

/* TODO once PeakIntensityRM are implemented 
	SpectrumLcms* spectrum=new SpectrumLcms(t);
	spectrum->setDescription(title);
	for (int i = 1; i <= sppts; ++i) {
		double pit = spec->GetDataPointYValue(i);
		
		if (! (pit > 0.)) continue; // skip if there's no peak intensity
		double pmz = 0., fmz = 0.;
		spec->GetQ1Mass(i, &pmz);
		spec->GetQ3Mass(i, &fmz);

		spectrum->push_back(new PeakIntensityRM(pmz, fmz, 0, pit));
	}

	// congratulation, we'd created a virtual MS^1 spectrum
	run.msCollection().push_back(spectrum);
*/
	// first cache all Q1 (=precusor) mass settings
	// the Q3 (=fragment) are simply ignored
	double *mz = new double [sppts];
	double *it = new double [sppts];
	int csize = 0;
	if ((! skip) || (checkmap(curpe, curxp, curcy)))
		for (int i = 1; i <= sppts; ++i) {
			double pit = 0., pmz = 0.;
			try {
				dhCheck( dhGetValue(L"%e", &pit, spec, L"GetDataPointYValue %d", i) );
				if (! (pit > 0.)) continue; // skip if there's no peak intensity
				dhCheck( dhCallMethod(spec, L"GetQ1Mass %d,%&le", i, &pmz) );
			} catch (string errstr) {
				cerr << "Error getting peaks" << endl
					<< endl << "details:" << endl << errstr << endl
					<< endl << "skipping spectrum" << endl;
				return;
			}
			
			mz[csize] = pmz;
			it[csize] = pit;
			++csize;
		}

	// Store into Captain Bol
	// *** EXPORT CAPTAINBOL
	SpectrumLcms* spectrum=new SpectrumLcms(t);
	spectrum->setDescription(title);

	spectrum->setRefid(mkrefid());
	
	if ((! skip) || (checkmap(curpe, curxp, curcy)))
		for (int i = 0; i < csize; ++i) {
			// skip flagged M/Z
			if (mz[i] <= 0.) continue;
			
			double	pmz = mz[i],
				pit = it[i];
			// sum all mass close to pmz
			for (int k=i+1; k < sppts; ++k) {
				if (fabs(pmz - mz[k]) < 0.01) {	// ridiculously close delta. FIXME maybe DeltaMass should work better here...
					mz[k] = 0.;
					pit += it[k];
				}
			}
			
			// add the virtual MS^1 peak to the spectum
			if (pit > 0.) // (only if there's peak intensity)
				spectrum->push_back(new PeakIntensity(pmz, 0, pit));
		}

	// congratulation, we'd created a virtual MS^1 spectrum
	run->msCollection().push_back(spectrum);

	delete mz; delete it;
}






void analyst::do_single_spectrum_ms1(const int sam, const int pe, const int xp, const float t, const double tic, const bool skip, bool docentroid) {
	// MS level 1 - usually the experiment #0 in an IDA
	// used to detect precursors
	cout << "\t\tGetting MS1 spectrum " << curcy << "\tat " << t << "\ttic:" << tic;

	bool centroided = false;
	int sppts = 0;
	int scans = 0;
	int stype = 0;
	string title("");

	try {
		char *ptr = NULL;
		dhCheck( dhCallMethod(spec, L"SetSpectrum %d,%d,%d,%e,%e", sam,pe,xp,t,t) ); // at the same exact time as the point on the chromatogram
		// note : whatever the default time scale in the wiff file, the "SetSpectrum" method *needs* seconds.
#if DEBUG_DUMP_CSV
		string debugout ("out_" + boost::lexical_cast<string> (curxp) + "_spec_" + boost::lexical_cast<string> (curcy) + ".csv");
		dhCallMethod(spec, L"SaveToTextFile %s", debugout.c_str());
#endif

		dhCheck( dhGetValue(L"%d", &sppts, spec, L"GetNumberOfDataPoints") );
		dhCheck( dhGetValue(L"%d", &scans, spec, L"GetNumberOfScans") );	// = Number of different Scans used to produce this "spectrumdata" object
		dhCheck( dhGetValue(L"%d", &centroided, spec, L"IsCentroided") );
		dhCheck( dhGetValue(L"%d", &stype, spec, L"SpectrumType") );
		dhCheck( dhGetValue(L"%s", &ptr, spec, L"DataTitle") );
		title = ptr;
		dhFreeString(ptr);
	} catch (string errstr) {
		cerr << "Error getting spectrum" << endl
				<< endl << "details:" << endl << errstr << endl
				<< endl << "skipping spectrum" << endl;
		return;
	}

	if (scans > 1) {
		cerr << "we should use data obtain from only 1 spectrum !" << endl;
	}

	if ((! centroided) && docentroid) {
		// FIXME do centroid
	}


	cout	<< "\tnum of points: " << sppts
//		<< "\tspectra:" << scans	
			// Well we just have to pay attention that no more than 1 scan are used in the current spectra
		 	// (unless combining several scan in one on purpose)
			// TODO optionally we could add a few asserts but it's not required
		<< "\t" << (centroided  ? "centroided" : "raw")
		<< endl;

	
	// *** EXPORT CAPTAINBOL
	SpectrumLcms* spectrum=new SpectrumLcms(t);
	spectrum->setDescription(title);

	spectrum->setRefid(mkrefid());
	
	if ((! skip) || (checkmap(curpe, curxp, curcy)))
		getdata(spectrum, sppts);

	/* // currently ranks aren't needed in a loader

	SpectrumPeakIntensityUtilities::computeIntensityRank(*spectrum);
	sort(spectrum->begin(), spectrum->end(), Peak::ltptr);

	// to get spectra numbers
	cout << run.msCollection().size() << endl;
	*/
	run->msCollection().push_back(spectrum);
}


void analyst::do_single_spectrum_ms2(const int sam, const int pe, const int xp, const float t, const double tic, const bool skip, bool docentroid) {
	// MS level 1 - usually the experiment #0 in an IDA
	// used to detect precursors
	cout << "\t\tGetting MS2 spectrum " << curcy << "\tat " << t  << "\ttic:" << tic;

	bool centroided = false;
	int sppts = 0;
	int scans = 0;
	int stype = 0;
	string title("");
	double	it = tic, // initialise intensity with some bogus data
 			// NOTE using the TIC value *IS NOT* conceptually  correct
  			// but still provides some remote approximate of variation of precursor over time
 			// (SmileMS has been using a similar approach to build a precursor elution curve)
 			// Nonetheless the *actual* precursor intensity is set in the next step in buildptr :
 			// once we have all spectra loaded (MS1 and MS2) a precursor intensity is just a matter 
 			// of copying the value from one spectra to the other.
		MoZ = 0.0, ce = 0.0; // ce = collision enery
	unsigned short chval = 0;	// chval = charge state, as a numeric value
	short type = SCAN_UNKNOWN;
	
	try {
		double m2 = 0.;

		char *ptr = NULL;
		dhCheck( dhCallMethod(spec, L"SetSpectrum %d,%d,%d,%e,%e", sam,pe,xp,t,t) ); // at the same exact time as the point on the chromatogram
		// note : whatever the default time scale in the wiff file, the "SetSpectrum" method *needs* seconds.
#if DEBUG_DUMP_CSV
		string debugout ("out_" + boost::lexical_cast<string> (curxp) + "_spec_" + boost::lexical_cast<string> (curcy) + ".csv");
		dhCallMethod(spec, L"SaveToTextFile %s", debugout.c_str());
#endif

		dhCheck( dhGetValue(L"%d", &sppts, spec, L"GetNumberOfDataPoints") );
		dhCheck( dhGetValue(L"%d", &scans, spec, L"GetNumberOfScans") );	// = Number of different Scans used to produce this "spectrumdata" object
		dhCheck( dhGetValue(L"%d", &centroided, spec, L"IsCentroided") );
		dhCheck( dhGetValue(L"%d", &stype, spec, L"SpectrumType") );

		// NOTE the following Method isn't documented in the CookBook, but was obtained by speaking with the technical support of Applied
		dhCheck( dhCallMethod(wiff, L"GetScanInfo %d,%d,%d,%d,%&hd,%&le,%&le,%&le,%&hd",
				sam,pe,xp,curcy,// Warning: cycles are 0 to N-1 (not 1 to N like DataPoints)
				&type,&MoZ,&m2,&ce,&chval) ); // at the same exact time as the point on the chromatogram

		dhCheck( dhGetValue(L"%s", &ptr, spec, L"DataTitle") );
		title = ptr;
		dhFreeString(ptr);
	} catch (string errstr) {
		cerr << "Error getting spectrum" << endl
				<< endl << "details:" << endl << errstr << endl
				<< endl << "skipping spectrum" << endl;
		return;
	}

	cout << "\tce:" << ce << " ch:" << chval ;
//	cout << " t:" << type; // Type isn't interesting we already know it from the outer loop.


	if (scans > 1) {
		cerr << "we should use data obtain from only 1 spectrum !" << endl;
	}

	if ((! centroided) && docentroid) {
		// FIXME do centroid
	}

	cout << "\tPre: " << MoZ
		<< "\tnum of points: " << sppts
//		<< "\tspectra:" << scans	
			// Well we just have to pay attention that no more than 1 scan are used in the current spectra
		 	// (unless combining several scan in one on purpose)
			// TODO optionally we could add a few asserts but it's not required
		<< "\t" << (centroided  ? "centroided" : "raw")
		<< endl;

	// *** EXPORT CAPTAINBOL
	// ch val -> ch mask
	const Pchargemask chmask = (0 == chval)
		? Peak::string2chargemask( "?" ) // TODO ask PAB if a default charge mask wouldn't be better than "?"
		: Peak::string2chargemask( lexical_cast<string>(chval).c_str() );

	SpectrumLcmsms* spectrum=new SpectrumLcmsms(MoZ, chmask, it);
	spectrum->setDescription(title);
	spectrum->setRetentionTime(t);
	spectrum->setCollisionEnergy(ce);

	spectrum->setRefid(mkrefid());

	if ((! skip) || (checkmap(curpe, curxp, curcy)))
		getdata(spectrum, sppts);

	run->msmsCollection().push_back(spectrum);
}


void analyst::getdata(Spectrum< PeakIntensity* > *sp, int sppts) {
	for (int i = 1; i <= sppts; ++i) {
		double pmz = 0., pit = 0.;
		try {
			dhCheck( dhCallMethod(spec, L"GetDataPoint %d,%&le,%&le", i, &pmz,&pit) );
		} catch (string errstr) {
			cerr << "Error getting peaks" << endl
				<< endl << "details:" << endl << errstr << endl
				<< endl << "skipping spectrum" << endl;
			return;
		}

		if (pit > 0.) // (only if there's peak intensity)

		// *** EXPORT CAPTAINBOL
		sp->push_back(new PeakIntensity(pmz, 0, pit));
	}
}

/**
 * populate references inside bol::runlcmsms
 * so that MS2 spectra have reference pointing to their respective MS1 parent and parent peak
 * and, If we have a parent peak, use its intensity as the precursor intensity of the MS2 spectra 
 * (instead of using TIC - which is wrong)
 */
void analyst::buildptr() {
	// *** EXPORT CAPTAINBOL
	cout << "building references..." << endl;
	for (SpectrumLcmsmsCollection::iterator i = run->msmsCollection().begin();
		i != run->msmsCollection().end();
		i++) {

		double it = 0.;
		SpectrumPeakRefDualPtr<PeakIntensity>* refptr = FindParentSpectrum(run->msCollection(), (*i), (*i)->moz(), &it);

		// mutate the intensity to the new value
		if (it > 0.)
			(*i)->setIntensity(it);

		if (refptr != NULL) {
			// we have a reference to a MS^1 spectrum & peak....
			SpectrumPeakRefPtr baseptr(refptr);
			(*i)->setSpectrumPeakRef(baseptr);	// put the actual precursor's intensity now that we know it
	
			// if everything is okay, now we can travel from the current spectrum to the parent intensity
			// (and both should be the same)
			if (refptr->peak()) // check if there's one
				assert((*i)->spectrumPeakRef()->peak()->intensity() == it);
		}
	}
	cout << endl;
}



/**
 * This function builds the maps necessary to find spectra by RefId
 *
 * @param detector current detector
 * @param ptr pointer to a vector to 64bits refid
 */
void analyst::initextract_map(ExtractList *ptr) {
	// first clear the old map if needed
	
	if (ptr) {
		RefMap = new PERIODMAP;
		
		for (ExtractList::iterator i = ptr->begin(); i < ptr->end(); i++) {
			// use the union to convert the 64bits refid into numbers we need
			REFID rf = { *i };
			
			if (RefMap->find(rf.th.period) == RefMap->end())
				(*RefMap)[rf.th.period] = new EXPMAP;
			
			if ((*RefMap)[rf.th.period]->find(rf.th.experiment) == (*RefMap)[rf.th.period]->end())
				(*((*RefMap)[rf.th.period]))[rf.th.experiment] = new CYCLEMAP;
			
			(*((*RefMap)[rf.th.period]))[rf.th.experiment]->insert(rf.th.cycle);
		}

#if 0 // debug code
		for (PERIODMAP::iterator i = RefMap->begin(); i != RefMap->end(); i++) {
			for (EXPMAP::iterator j = i->second->begin(); j != i->second->end(); j++) {
				for (CYCLEMAP::iterator k = j->second->begin(); k != j->second->end(); k++) {
					cout << "[" << i->first << "][" << j->first << "][" << *k << "]" << endl;
				}
			}
		}
#endif
	}
}


/**
 * This function checks if Period, Experiment (and Cycle) are in the list of elements to be processed
 * @param pe 
 * @param xp 
 * @param cy
 */
bool analyst::checkmap(int pe, int xp) {
	if (! RefMap)
		return false;
	
	if (RefMap->find(pe) == RefMap->end())
		return false;
			
	if ((*RefMap)[pe]->find(xp) == (*RefMap)[pe]->end())
		return false;
	
	return true;
}
bool analyst::checkmap(int pe, int xp, int cy) {
	if (! checkmap(pe, xp) )
		return false;
	
	if ( (*((*RefMap)[pe]))[xp]->find(cy) == (*((*RefMap)[pe]))[xp]->end())
		return false;
	
	return true;
}

}
}

